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Abstract 

By means of a Quadratic Diffusion Monte Carlo method we have performed a 
comparative analysis between the Aziz potential and a revised version of it. The results 
demonstrate that the new potential produces a better description of the equation 
of state for liquid '*He. In spite of the improvement in the description of derivative 
magnitudes of the energy, as the pressure or the compressibility, the energy per particle 
which comes from this new potential is lower than the experimental one. The inclusion 
of three-body interactions, which give a repulsive contribution to the potential energy, 
makes it feasible that the calculated energy comes close to the experimental result. 
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I Introduction 

Many-body techniques have achieved a high level of accuracy in the description of atomic 
^He and "^He, which constitute the most characteristic examples of quantum liquids. The 
theoretical approaches to the many-body problem can be classified in two large blocks de- 
pending on the use or non-use of stochastic procedures. Among the non-stochastic methods 
it is the variational framework Q combined with integral equations such as HNC, which 
has provided the best results in the knowledge of the ground state. Also, perturbation 
schemes constructed on correlated basis (correlated basis function theory -CBF- Q) has 
led one to extend this study to the lowest excited states. On the other hand, Monte Carlo 
(MC) methods Q give exact information, within some statistical uncertainities, on the 
ground state of bosonic systems both at zero and finite temperature. The initial constraint 
imposed by the use of a finite number of particles in MC simulations does not influence 
appreciably to the energetic properties. However, the structure properties at r ^ oo 
{k 0) related to long-range correlations are out of scope. 

The high agreement between the theoretical results and the experimental data is also 
linked to the well known interatomic interaction for He atoms (pair wise additive form). 
For the last ten years, the HFDHE2 potential proposed by Aziz et al. [Q] has allowed 
for reproducing the energetic and structure properties of liquid He quite well both in 
homogeneous and inhomogeneous phases fo], |^ ^. Despite of the accuracy of this 
pair-potential a renewed version of it (HFD-B(HE)) was published by Aziz et al. in 1987 

The new Aziz potential (hereafter referred to as Aziz II potential) was brought about 
as a consequence of several new theoretical and experimental results which appeared in 



the literature between the publication of the two potentials. First, Ceperley et al. [10| 



pointed out by means of a quantum Monte Carlo calculation of the interaction energy of 
two He atoms, with internuclear separations less than 1.8 A, that the Aziz potential is 
too repulsive below this distance. On the other hand, new experimental measurements of 
the second virial coefficients and transport properties for '^He and ^He showed evidence of 
some small inconsistences of the Aziz potential. The explicit expressions of the Aziz and 
Aziz II potentials appear in the Appendix A. Apart from a soft core, the Aziz II potential 
has its minimum at e = 10.95 K, = 2.963 A while Aziz potential has its minimum at 
e = 10.80 K, Vm = 2.967 A. Therefore, the new potential is only slightly deeper with the 
minimum localized at a lower interatomic separation. 

To start on a theoretical comparative study between He potentials it is necessary to 
calculate the properties of the liquid as precisely as it is possible. Stochastic methods 
provide the appropriate tools for this purpose, especially in the case of bosonic systems as 
^He. In the past, the Green's function Monte Carlo method (GFMC) was used to elucidate 
between different models for the pair interaction. The main conclusion of this analysis Q 
stated that the Aziz potential was the best interaction to study the properties of liquid 
and solid helium. 

Our objective in the present work is to perform a comparative analysis between the 
two Aziz potentials to establish if the new potential (Aziz II) produces even better results 
than the previous one. The calculation presented here follows an alternative procedure to 
GFMC known as Diffusion Monte Carlo (DMC). 

Both GFMC, developed by Kalos et al. (l|, ||], and DMC algorithms |l|, ||] solve 
stochastically the Schrodinger equation in imaginary time. The GFMC scheme constructs 
a time integrated Green's function by means of a double Monte Carlo sampling. On the 
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other hand, the DMC algorithm is a simpler method that assumes an approximate form 
for the Green's function for small time-steps At. In this case, after an iterative process and 
sufficiently long times, only the ground state wave function survives. Therefore, the exact 
energy per particle of the system is obtained when the limit At ^ is considered. DMC 
is posterior to GFMC but up to now it has already been applied to the study of small 
molecules |14|, solid hydrogen |15| or ^He clusters ||^]. The main disadvantage of the DMC 
algorithms used in the major part of those works is that the energy eigenvalues change 
linearly with At. This fact obliges one to perform several calculations using different 
values for the time-step and next to extrapolate the exact value in the limit At — s- 0. To 
avoid this difficulty several quadratic algorithms have been devised but the success of this 
improvement has not been complete. Recently, a new Quadratic Diffusion Monte Carlo 
(QDMC) method has proved to work efficiently in the description of ^He droplets 
0. In the present work, we use a QDMC method with a very similar algorithm to the 



one reported in Ref.|l(:]. In the next sections of the article we will justify the accuracy 
of the proposed method which allows for the possibility of calculating the properties of 
the system at a finite time-step without introducing any significant difference with the 
extrapolated value. 

The outline of this paper is as follows: In Sec. II the Quadratic Diffusion Monte Carlo 
method to solve the Schrodinger equation is presented. The consistency of the algorithm 
is checked by using different trial functions and several numbers of particles. The time- 
step dependence of the energy per particle shows the expected quadratic behaviour. A 
comparative analysis of the two Aziz potentials is reported in Sec. III. A perturbative 
estimation of the contributions coming from various three-body potentials is also reported. 
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A brief discussion and conclusions comprise Sec. IV. 



II Computational algorithm 

The starting point in Diffusion Monte Carlo methods is the Schrodinger equation for 
particles written in imaginary time: 



-^^ = {H-E)^{R,t) (1) 



where R = (ri, . . . , fat) and t is measured in units of h. ^'(R, t) can be expanded in terms 
of a complete set of eigenfunctions (/>i(R) of the Hamiltonian: 

^(R,t)=Y,Cn exp[-{Ei-E)t] <j)i{-R) , (2) 

n 

where Ei is the eigenvalue associated to 0j(R). The asymptotic solution of Eq. (|^) for any 
value E close to the energy of the ground state and for long times {t — > oo) gives 00 (R), 
provided that there is a non-zero overlap between ^'(R, t = 0) and the ground state wave 
function i;^o(R-)- 

In a computer simulation of Eq. (|l|) it is crucial to use the importance sampling tech- 



nique |12| in order to reduce the statistical fluctuations to a manageable level. Following 



this method, one rewrites the Schrodinger equation for the function: 

/(R,t) =V(R)^(R,t) , (3) 
where ip(R-) is a time-independent trial function. Considering a Hamiltonian of the form: 
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Eq. (0) turns out to be: 



^ffli^ = -DVyi-R,t)+DVu{F{'R)f{'R,t)) + {ELi-R)-E) f{-R,t) 



= {Ai + A2 + A3)f{K,t)=Af{K,t), (5) 
where D = 'h? /{2m), E^iYC) = '4){YC)~^Hil:{R) is the local energy, and 

F(R) =2V(R)-^VrV^(R) (6) 

is called the drift force. F(R) acts as an external force which guides the diffusion process, 
involved by the first term in Eq. (^), to regions where V'(r^) is large. 
The formal solution of Eq. (^) is 

/(R', t + M)= J G(R', R, At) /(R, t) dR (7) 

with 

G(R',R,At) = (R'l exp(-^At) [R) . (8) 

While GFMC method works with the whole Green's function, DMC algorithms rely 
on reasonable approximations of G(R', R, At) for small values of the time-step At. Then, 
Eq. (^) is not directly solved but iterated repeatedly to obtain the asymptotic solution 
/(R,t^cx)). 

In the Quadratic Diffusion Monte Carlo algorithm we have used, the Green's function 
G(R',R, At) is approximated by: 



exp {-A At) 



(9) 



exp ^-^3^^ exp ^-^2^^ exp(-y4iAt) exp ^-^2^^ exp (^-A^ 



At 



This decomposition, which is not unique [16|, is exact up to order (At)^. Assuming (9), 
Eq. (0) becomes: 



f(K',t + At) = J 



with 



G3 (R',Ri,^)g2 (^Ri,R2,^) Gi(R2,R3,At) (10) 



G2 [ R3, R4, — G3 ^R4, R, — 



/(R,t)dRi ...dR4dR, 



Gi (R', R, t) 



3N 

{A-TTDt) 2 exp 



fR' - R) 



ADt 



and 



R(0) = R 



G2(R',R,t) = 6{R'-K{t)), where < 



^ d^ = DF{K{t)), 



GsiK', R, t) = exp [-(^l(R) - E) t] 6{R' - R). 



(12) 



(13) 



In our Monte Carlo computations, /(R, t) is represented by walkers Rj , each one 
representing a set of the 3A^ coordinates of the particles. The algorithm used for the 
implementation of Eq.(lO) goes through the following steps: 



1. Move the walkers, under the drift force F(R), during an interval At/2 with accuracy 
(At)2. 



2. Apply to each walker a displacement x randomly drawn from the 3iV Gaussian 
distribution exp(-x^/(4DAt)). 

3. Repeat step 1. 

4. Randomly replicate each walker rir times, in such a way that 



{rir) = exp 



5. Go to step 1 for the next walker Rj , until the set of walkers is exhausted. The new 
set obtained corresponds to /(R, t + At). 

The whole procedure is repeated as many times as it is needed to reach the asymptotic 
limit {t oo). From then on, the walkers Rj are used to obtain the expectation values of 
the magnitudes to be determined. 

In order to establish the preciseness of the method several aspects have to be consid- 
ered. First, Monte Carlo information about /(R, t) only allows measurements of quantities 
by means of mixed estimators, i.e., Thereby, only when the operator A is the 

Hamiltonian the mixed estimator gives the exact expectation value. To obtain other 



ground state properties a simple linear extrapolation |11] has been currently used: 



(^'1^1^') = 2(V^|A|^') - (V'l^lV') (14) 

This method involves the performance of a Variational Monte Carlo (VMC) calculation 
to determine the variational expectation value (iplAlip) . It is interesting to notice that 
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a VMC calculation can be carried out with the same algorithm described for DMC only 
suppressing the branching term (p^). 

In Table I, the results for the potential, kinetic and total energies per particle obtained 
with VMC and QDMC methods are shown together with a GFMC result They corre- 
spond to an Aziz potential calculation with N = 128 at density p = 0.365 a^^ {a = 2.556 
A). The trial wave functions ^/^ji and contain different two-body correlation factors, 
and ipjT includes also three-body correlations. Explicit expressions of these trial functions, 
together with the values of the parameters involved in them, are given in Appendix B. 

As it is shown in Table I, there are not significative discrepancies between the QDMC 
results for the total energy. The perfect agreement between the QDMC results and the 
GFMC value is also remarkable. Equation (^) is used to estimate the kinetic and potential 
contributions to the total energy. In spite ot its simplicity, this method gives very similar 
values for the partial energies even when trial wave functions as different as the ones 
reported in Table I are used as importance sampling. New methods to avoid the slight 
influence of the trial wave function in the extrapolated estimators have been recently 



suggested by Barnett et al. [T^] and Zhang and Kalos |18]. 

The effect of a finite volume simulation box has also been considered, raising the 
number of particles N from N = 128 (which has been used for the bulk of the calculation) 
up to = 190 for and p = 0.365 a^^. The differences encountered were compatible 
with the size of the statistical fluctuations reported in Table I. 

Another important parameter in the calculation is the population of walkers nw All 
the results reported in the present work have been obtained after a preliminary analysis of 
the influence of the population in the average energy of the system, and the final results 
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correpond to the empirical asymptotic value of ra^. This asymptotic population decreases 
when one improves the quality of the trial wave function ijj. In fact, whereas Tiyj — 400 
for ipji and ipj2 this value is reduced to = 250 in the ipjT case. Actually, we have 
developed a parallel QDMC code based on the equal role played by the different walkers. 
The calculation has been carried out on a massively parallel computer CM2, which best 
performance is obtained when a large number of walkers is considered. The length of the 
series have been 12000-15000 for ^ji and and 8000-10000 for Vjt- 

A final but not less important point is the time-step dependence of the QDMC algo- 
rithm. In Fig. 1 it is shown a characteristic result of the total energy as a function of the 
time-step At. The time t is measured in reduced units r, where 



As one can see in the Figure, there is a clear departure from the linear time dependence 
supplied by the linear DMC algorithms. If a second order polynomial fit E/N = {E/N)q + 
A(At)^ (solid line in the Figure) to the QDMC results is performed, one obtains an 
extrapolated value of {E/N)q = (—7.124 zt 0.003) K, which is indistinguishable from the 
values obtained working with At = (1 — 2) ■ 10~^ r. Therefore, it is plausible to calculate 
the properties of the system accurately using a single value for At, lying in the stated 
range, without the necessity of a complete analysis in time to extrapolate the correct 
results. 



9 



Ill Results 

In this section the numerical results for the energy and for the structure properties using 
the Aziz and Aziz II potentials are presented. First, we analyse the differences between 
the two interatomic potentials and then, the contribution to the total energy of several 
models for the three-body interactions. In all the calculations reported below the trial 
function i{jj2 (see Appendix B) has been used as importance sampling. The average 
population size (riw) ranges from 400, near the equilibrium density {p = 0.365 a~^), to 
900 for the highest densities. On the other hand, the value for the time-step has been 
taken as At = 1.25 • 10^'^ r around the equilibrium density and At = 1.0 • 10^^ r for higher 
densities. No significant deviations in the results of the energy are observed when these 
At values are doubled. 

A. Two-body potentials: Aziz vs. Aziz II 

As it has been commented in the Introduction, the differences between the Aziz and 
Aziz II potentials are not very large. However, slight differences in the values of the 
parameters entering into V{r) produce relatively large changes in the energy as it was 
asserted by Kalos et al. ||5| in their GFMC calculation of the equation of state of liquid 
■^He using the Aziz potential. The energies obtained for both potentials, together with 



the experimental results of Ref. [19| are reported in Table II. In parenthesis there are 
the GFMC results for the Aziz potential. The GFMC and QDMC calculations are in 
good agreement, but a small deviation between both results is observed at high densities. 
The kinetic and potential energies are also given in the Table. The potential energy has 
been calculated by means of the extrapolated estimator (Eq. 0), and the kinetic energy 
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comes from the difference between the total and the potential energies. A comparison 
between the partial energies of the two potentials reveals that, while the kinetic energy is 
practically the same, the Aziz II potential energy is, in absolute value, larger than the Aziz 
one. In particular, the Aziz II potential lowers the potential energy with respect to the 
Aziz case in a quantity which grows from ~ 0.19 K at p = 0.365 to ~ 0.23 K at the 
highest density p = 0.490 a^^. The partial energies for both potentials satisfy the lower 
bound for the kinetic energy and the upper bound for the potential energy (T/N > 13.4 K 
and V/N < —20.6 K, at the equilibrium density) [pO| . 

Concerning the total energies, one can observe that the experimental values are ap- 
proximately located at the middle of the Aziz and Aziz II results. This fact is clear from 
Fig. 2 where the equation of state of liquid ^He is shown in comparison with the experi- 
mental results. The lines in the Figure correspond to numerical fits to the results reported 
in Table II, excluding the highest density (0.490 cr~^) because it is quite far from the 
experimental freezing density pi = 0.430 a~^. In the majority of microscopic calculations 
on liquid He it has been used a polynomial fit of the form 



where e = E/N and po is the equilibrium density, to determine the equation of state. On 
the other hand, in calculations based on Density Functional Theory the form 




(16) 



e = h p + cp 



1+7 



(17) 



proposed by Stringari and Treiner |21], has proved to be very efficient in describing prop- 
erties of homogeneous and inhomogeneous (including an additional surface term in Eq. 
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^) liquid ^He. In the case of the Aziz II potential both analytic forms are compatible 
with the results of the energy, taking into account their respective errors. However, we 
have checked that only the function (^) provides the correct result for the energy when 
densities lower than 0.328 are considered. Therefore, all the results presented below, 
concerning the equation of state, are derived starting on the second option (Eq. p!7|) . 
The values of the parameters which best fit our Aziz II potential results are: 

b = (-27.258 ±0.017) Ka^ 
c = (114.95 ± 0.22) 

7 = 2.7324 ±0.0020. (18) 

The same analysis has been performed by taking the energy results of the Aziz poten- 
tial. In this case, neither the polynomial form (|l6| ) nor the Stringari's ( |T7|) are statistically 
compatible with our results. This fact is clearly reflected in Fig. 2, where several Aziz 
points (the size of each point is larger than its error bar) are not intersected by the result 
of the fit (represented with a dashed line). In spite of this severe restriction, and to make 
possible the comparison with the equation of state provided by the Aziz II potential, the 
optimum values 

b = (-26.947 ± 0.016) Ka^ 
c = (115.72 ±0.21) Kcr3(^+^) 

7 = 2.7160 ± 0.0020 (19) 

are taken. 
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We have also fitted the same type of function to the experimental results of Ref. [19|. 
In this case, the parameters b and c have been fixed to reproduce the equilibrium density 
and the energy at this density, whereas the parameter 7 has been obtained by means of a 
numerical fit to all the energies reported in that work. The values obtained are 

b = -26.746 Ka^ 

c = 116.69 E:o-^(^+t) 



7 = 2.7773. (20) 

Once the equation of state e{p) is known, it is straightforward to calculate the isother- 
mal compressibility, defined as 

where P{p) = {de /dp) is the pressure, and the velocity of sound given by 

c{p)={ . (22) 

In Table III the results of the pressure, the compressibility and the velocity of sound of the 
two Aziz potentials are compared with the corresponding experimental values at the exper- 
imental equilibrium density (/Og^^ = 0.3646 o""'^). It is remarkable the accuracy provided 
by the Aziz II potential, giving results for these quantities which are indistinguishable 
from the experiment. Conversely, the equation of state corresponding to the Aziz poten- 
tial supplies results which are slightly worse. The differences between the equations of 
state for the two potentials remain when the density increases, as one can see for P{p) in 
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Fig. 3 and for k{p) in Fig. 4. The equation of state corresponding to the Aziz II potential 
gives an excellent description of these magnitudes for all the values of the density here 
considered. 

Apart from the ground state energy, the Monte Carlo methods yield other interesting 
information. The radial distribution function 

N{N-1) /|^(ri,...,rAf)|^ dr^-.-drN 
'^'''^ = /l^(r„...,r.)|^ dr,...dr, ^^'^ 

and its Fourier transform, the static structure function 

S{k) = l + p J dr e^^ "" {g{r) - 1) (24) 

are fundamental in the study of fluids. The calculation of these quantities is more involved 
than the calculation of the energy [ pT| , p^ ], but the extrapolation procedure (Eq. |l^ allows 
results which are practically independent of the trial function used as importance sampling. 

The radial distribution function g{r), obtained in a Aziz II calculation at a density 
p = 0.365 (7^^, is shown in Fig. 5 in comparison with an experimental determination at 



T = 1.0 K hy Svensson et al. [23|. There is a good agreement between the calculated and 
the experimental g{r), mainly in the first peak. In Fig. 6, the structure function S{k)^ 
obtained by means of a Fourier transform of the g{r) shown in Fig. 5, is plotted together 
with the experimental measure of Ref. ||2^. Due to the finite size of the simulation box, 
there are not reliable results for S{k,) for A; ^ 1 A^^. The theoretical S{k) is again very 
close to the experimental result, but the height of the experimental main peak is slightly 
higher. On the other hand, other experimental determinations of S{k) |^4| point to lower 
values of the intensity of the first peak, even below our results. In fact, analysis of the 
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influence of the temperature T in S{k) [p4| , 25 1 indicate that the largest variation of the 
structure function with T is placed in the vicinity of the flrst peak. 
The one-body density matrix p{r) deflned as 

/ ^(r^i , . . . , rjv) ^(ri , ■ ■ ■ , rjv) dra ■ . . drN , . 

piriv) = N ^— — r|2— ^ ^ 25 

and its Fourier transform, the momentum distribution 

n{k) = (27r)3 puQ 5{k) + p j dr e^^ '' {p{r) - p{oo)) (26) 

can also be computed using the conflgurations generated by the QDMC code. The function 
p(r) is obtained as the expectation value of the operator 

*(ri,---,r, + r,...,r^) ^ ^27) 



^(ri,...,r7v) 

evaluated on the configuration space over a set of random desplacements of the particle i. 
The condensate fraction no, i.e., the fraction of particles occupying the zero momentum 
state, may be extracted from p{r) by means of the asymptotic condition 

no = lim p{r) . (28) 



In Fig. 7 the momentum distribution obtained via the Eq. (26) is plotted, as kn{k), 
for three values of the density. The correlations between the particles make the population 
of states with high momenta increase with the density. The shoulder observed at A; ~ 2 
A^^ for the three curves, which has been observed in other theoretical calculations of 



n( 



(k) |22, 27|, has been attributed in the past to the zero-point motion of the rotons [p8| . 
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On the other hand, it has been proved that if the condensate fraction is non-zero, n{k) 
diverges as 1 / A; when A; ^ [^]. Again, the finite value of the simulation cell precludes 
the possibility of reproducing this behaviour. 

We have also determined the condensate fraction from the extrapolated estimation of 



p{r) and the relation (25). At the equilibrium density, we get no = 0.084 it 0.001 which 
is a value slightly smaller than the one obtained in a GFMC calculation (0.092 it 0.001) 
[ p^ using the Aziz potential. The discrepancy between the two results are not due to the 
use of different potentials. In fact, we have calculated p{r) for the two Aziz potentials and 
no significant differences appear. The same conclusion holds for the radial distribution 
function g{r). 

A final point of interest is the density dependence of the condensate fraction . In Fig. 
8, the change in the value of uq is shown for a wide range of densities. The condensate 
fraction decreases with the density, following a law nearly quadratic in p. In the Figure, 
a quadratic fit to the results is shown as a "guide to the eye" . 

B. Three-body interactions 

The importance of three-body interactions in helium has been discussed for a long 
time. It has been argued that these interactions would be present in He but its relative 
contribution to the total energy is still open to question. The most widely known model 
for the three-body potential is the triple-dipole interaction derived by Axilrod and Teller 
considering perturbative theory. The Axilrod- Teller (AT) potential, which has been 
usually considered as the major contribution to the energy coming from the three-body 
interactions, provides a positive correction to the potential energy. The amount of this 
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effect was calculated for the first time by Murphy and Barker |31] by means of a Vari- 
ational Monte Carlo calculation. Afterwards, that contribution was estimated by Kalos 
in a Rayleigh-Schrodinger perturbative calculation starting on GFMC con- 
figurations. From this analysis it was pointed out that, on the one hand, the GFMC 
prediction for the expectation value of the three-body potential (V3) was in accordance 
with the Murphy and Barker's variational results and, on the other, there were no relevant 
differences between the results coming from a Lennard-Jones and an Aziz potential cal- 
culations. Another conclusion of these GFMC works was that the inclusion of three-body 
potential contributions on the total energy worsened the two-body results along the whole 
equation of state. 

In spite of the AT potential being the dominant contribution to (V3), it has been 
proved that at short interparticle separations a non-additive and attractive force emerges. 
This short-ranged three-body interaction , usually known as exchange interaction, is due 
to the influence in the charge densities of two interacting atoms by the presence of a third 



near particle. Bruch and McGee [32| proposed a model potential (BM) to account for this 
effect, fitting the parameters of the exchange part to the atomic calculations of the energy 
of three He atoms at very short distances from Novaro and Beltran-Lopez |3^. Loubeyre 

has proved that the BM three-body potential, in conjunction with the Aziz potential, 
accurately describes solid helium at high pressures and room temperature. The explicit 
forms for the AT and BM potentials are given in Appendix A. 

As it has been previously discussed, the Aziz II results for the energy per particle are 
below the experimental results for all the densities considered (see Fig. 2). Therefore, 
the inclusion of a repulsive contribution to the potential energy, arising from three-body 
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interactions, could bring the theoretical results nearer to the experimental ones. In Table 
IV the results for the total {E/N) and potential {V/N) energies are reported in compar- 
ison with the experimental values of the energy. In all cases, the three-body potential 
energy is obtained by means of a Rayleigh-Schrodinger perturbative calculation, following 



the method described by Kalos et al. [22]. As one can see, the AT potential produces an 
increase in the energy, leading to values which are slightly higher than the experiment. 
Moreover, the difference between the Aziz II-I-AT and experimental values increases ap- 
preciably with the density, yielding to poor results for derivative magnitudes of the energy 
as the pressure or the compressibility. The results of the energy, using the BM poten- 
tial, appear in the second column of Table IV. The exchange part of the BM potential 
practically cancells the repulsive contribution of the dispersion term (AT) becoming even 
dominant at the highest densities. The resulting energies lie very near to the two-body 
calculation but also in this case, as in the AT one, with a worsening reproduction of the 
dependence of the energy with the density. Therefore, neither the simple AT potential nor 
the more elaborated one (BM) improve, in a significant way, the Aziz II results. In fact, 
it seems more convincing that, in the density regime of liquid ^He, the main three-body 
contribution comes from the AT potential, the exchange part of the BM potential being 
too large. We should notice that the parameters of the BM potential have been fitted 
to reproduce the energy of helium trimers with interparticle separations considerably less 
than the characteristic mean distance between the atoms in the liquid. Then, it is uncer- 
tain that the same parameters, or even the same analytical form, could be used to study 
the liquid phase. 
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In the third column of Table IV, labelled as MBM, we report the results which are 
obtained by using the BM potential with a modified value A' = A / 3 (see Appendix A) . 
Now, the energy at the experimental equilibrium density reproduces the experimental 
result and a quite good description is also obtained at higher densities. In Fig. 9, the 
equation of state obtained with the Aziz II+MBM model is depicted together with the 
experimental and the Aziz II results. The values of the parameters of the fit for the Aziz 
II+MBM calculation are: 

b = (-27.202 ± 0.017) Ka^ 
c = (114.11 ± 0.21) 

7 = 2.6961 ± 0.0020 (29) 

The Aziz II+MBM results for the pressure and the compressibility are plotted in Fig. 
10 and Fig. 11, respectively, in comparison with the experimental values. One can observe 
that there are slight differences between the theoretical and the experimental results, which 
are more evident in the pressure case. In fact, these discrepancies reflect the departure of 
the Aziz II+MBM total energies from the experimental values when the density increases. 
This small effect on the energy, which can be observed in Fig. 9, is enlarged when the 
derivative magnitudes of the energy as P{p) or /c(p) are calculated. 

IV Discussion 

The properties of bulk liquid ^He have been investigated by means of the Diffusion Monte 
Carlo (DMC) method. It has been proved that the extension of the DMC algorithm up 
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to second order (QDMC) allows for the possibility of calculating the energy without the 
extrapolation to At = 0, required in the linear DMC codes. We have applied the QDMC 
method in order to perform a comparative analysis between the Aziz and the Aziz II 
two-body potentials. The calculations have been extended to a wide range of densities 
in order to contrast the theoretical predictions on the equation of state provided by the 
two Aziz potentials. The results unambiguously demonstrate that the new Aziz potential 
gives better results than the old one, especially when the dependence of the pressure and 
the compressibility on the density is considered. In particular, the Aziz II results for P{p) 
and k{p) are indistinguishable from the experimental values. However, the results for 
the energy are below the experimental determinations. This difference could suggest the 
presence of three-body interactions in He. 

We have performed a Rayleigh-Schrodinger perturbative estimation of the three-body 
potential energy using two different models. The results obtained have shown that nei- 
ther the triple-dipole potential of Axilrod- Teller nor the Brunch-McGee potential, which 
includes the exchange interaction at short distances, improve the equation of state given 
by the Aziz II potential. To make the three-body correction compatible with the experi- 
mental results a simple change in the parameters entering into the BM potential has been 
examined (MBM potential) . The Aziz II-I-MBM model provides a good description of the 
equation of state {E / N){p) but the results for P{p) and worsen with respect to the 
ones calculated with the Aziz II only. On the other hand, we would point out that the 
Aziz II results for the energy are shifted with respect to the experiment in a constant value 
for all the densities. This fact explains the excellent description of P{p) and k{p) given by 
the new Aziz potential. The conclusion is that to account for the experimental energies. 
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a constant value for (V3) would be required, although from the theoretical point of view 
it seems more plausible a correction which becomes larger when the density increases. 

Concerning other properties as the radial distribution function or the momentum dis- 
tribution, no significant differences between the results given by the two Aziz potentials 
are observed. Overall, the agreement between the QDMC results and the experiment is 
quite satisfactory. Finally, we would remark that the accuracy of the Aziz II potential 
in describing the bulk ^He liquid phase makes it recommendable for future calculations 
of the solid phase or films, especially when the derivative magnitudes of the energy are 
among the main objectives. 
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Appendix A: Two- and three-body potentials 



The form of the HFDHE2 (Aziz) potential ||] is 



V{r) = e 



A exp(-ax) - F{x) ^ {C2j+& / x 

j=0 



2j+e 



where 



F(x) 



exp 



X < D 
x> D 



with 



The values of the parameters for the Aziz potential are 



e = 10.8 K Cg = 1.3732412 



= 2.9673 A Cs = 0.4253785 



D = 1.241314 Cio = 0.1781 



a = 13.353384 A = 0.5448504 • 10^ 



(A.l) 



(A.2) 



(A.3) 



(A.4) 



The HFD-B(HE) (Aziz II) potential Q, which is quite similar in form to the Aziz 
potential, is given by 



V{r) = e 



Z 

A exp (-ax + /Jx^) - F(x) ^ {C2j+6 / '■ 



.2j+6 



j=0 



(A.5) 



where the function F(x) and x are formally the same as in the Aziz potential (Eqs. 



A. 2, A. 3). The values of the parameters for the Aziz II potential are 
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£ = 10.948 K 



Ce = 1.36745214 



r™ = 2.963 A 



Ch = 0.42123807 



D = 1.4826 



Cio = 0.17473318 



(A.6) 



a = 10.43329537 A = 1.8443101 • 10^ 
(5 = -2.27965105. 

The models for the three-body interactions we have used are those given by the Axilrod- 
Teher (AT) @ and Brunch-McGee (BM) ||] potentials. The form of the AT potential 
is 



V3(r-i2,ri3,r23) 



1/(1 + 3 cos 4>i cos 4>2 cos I 

3 3 3 

fy-t^-> iyi^J ly^^ 

'12 '13 '23 



(A.7) 



where (j)2 and (/>3 are the interior angles of the triangle formed by the three atoms. We 
use the Leonard's helium value v = 0.327 Ka'^ [^], assuming the radial distances rij in a 
unities. 

The BM potential is given by 



V3(n2,n3,?^23) 



'"is ^"23 



A exp (-a (ri2 + ris + r23) 



X (1 + 3 cos </>i cos 4>2 COS 03 ) 



(A. 



The values for the two new parameters appearing in the BM potential are 



A = 9676545.53 K a = 4.948 a'^ . 



(A.9) 
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Appendix B: Trial functions 



In this Appendix we give the exphcit forms of the trial functions used as importance 
samphng in the QDMC calculations as well as the values of the parameters involved. The 



first one is the well known McMillan two-body trial function |36 



V'ji = n 

i<j 



11 

2 Ir 



(B.l) 



We have taken the value b = 1.20 a which optimizes the VMC energy at the experimental 
equilibrium density. 

Most of the present work has been carried out using the Reatto two-body function [37| 



V' J2 = n 

i<j 



1 / b 



L 



'2 UiJ 2 ^""P 



A 



(B.2) 



with L = 0.2, A = 2.0 cr, A = 0.6 a and b = 1.20 a. These values, optimal at the 
experimental equilibrium density, have also been used for the other densities. 

The third trial function, which was proposed by Schmidt et al. [^, contains two- and 
three-body correlations. It is explicitely given by 



1 



1 



where 



l^k 



and 



i{r) = exp 



r -n 



(B.3) 



(B.4) 



(B.5) 
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The values for the triplet parameters, roughly optimal at the equilibrium density, are: 
A = -1.08 n = 0.82 a and = 0.50 a. 
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Table Captions 

Table I: Results for the total, kinetic and potential energies for different trial wave func- 
tions. The forms of ipji, ipj2 and ^jt as well as the values of the parameters entering into 
them are noted explicitely in Appendix B. In the last row, the GFMC results from Ref. 
[H (a) and Ref. (b) are also reported. All energies are in kelvin per particle. 

Table II: Results for the total and partial energies from the QDMC calculations with 
the Aziz potential, the Aziz II potential and experiment [|^. The numbers quoted in 
parenthesis are taken from Ref. |^]. All energies are in kelvin per particle. 

Table III: QDMC results for the pressure P, the compressibility k and the velocity of 
sound c at the experimental equilibrium density using the Aziz and Aziz II potentials. 
The last row contains the experimental values derived from the experimental equation of 
state (20). 

Table IV: Energies from the QDMC calculations with the Aziz II potential including 
the perturbative estimation of the expectation value of several models for the three-body 
interactions (AT, BM, MBM; see text). The last column contains the experimental values. 
All energies are in kelvin per particle. 
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Figure Captions 

Fig. 1. Time-step dependence in the QDMC method. The solid Hne is a second order 
polynomial fit to the calculated points. 

Fig. 2. Equation of state for liquid ^He. The circles are the QDMC results with the 
Aziz potential and the dashed line is a fit to the calculated energies. The solid circles 
correspond to the QDMC energies with the Aziz II potential; the solid line is a fit to 
these energies. These fits have been performed with Eq. (17), the parameters being those 
given in (18) and (19) for the Aziz II and Aziz potentials, respectively. The experimental 
values, represented by solid triangles, are taken from Ref. [^]. The error bars of the 
QDMC results are smaller than the size of the symbols. 

Fig. 3. Pressure of liquid ''He as a function of the density. The circles and solid circles 
correspond to QDMC calculations with the Aziz and Aziz II potentials, respectively. The 
dashed and solid lines are numerical fits to the Aziz and Aziz II pressures, respectively. 
The experimental results, from the experimental equation of state (20), are represented 
by diamonds which are practically hidden below the Aziz II values. The error bars of the 
QDMC results are smaller than the size of the symbols. 

Fig. 4. Isothermal compressibility of liquid ^He as a function of the density. The same 
notation as in Fig. 3. 

Fig. 5. Two-body radial distribution function at the experimental equilibrium density. 
The solid line is the QDMC result and the solid circles correspond to the neutron diffraction 



experimental determination from Ref. [23|. 



30 



Fig. 6. Static structure function at the experimental equilibrium density. The solid line 
is the QDMC result, obtained by a Fourier transform of the radial distribution function 
showed in Fig. 5. The solid circles are the experimental determination from Ref. [23|. 



Fig. 7. Dependence of the calculated momentum distribution on density. The long- 
dashed, solid and short-dashed lines stand for the results at densities of 0.328 cr~^, 0.365 
and 0.401 o"~^, respectively. 

Fig. 8. Condensate fraction in liquid ^He as a function of density. The solid line is a 
second order polynomial fit to the calculated values. The error bars of the results are 
smaller than the size of the symbols. 

Fig. 9. Equation of state of liquid ^He. The circles are the QDMC results with the 
Aziz II+MBM potentials; the solid line is a fit to these energies. The solid circles and the 
dashed line correspond to the calculation with the Aziz II potential. The solid triangles 
are the experimental values form Ref. |jl^. The error bars are smaller than the size of the 
symbols. 

Fig. 10. Pressure of liquid ^He as a function of the density. The solid circles are the Aziz 
II+MBM results and the diamonds are the experimental values from the equation of state 
(20). The solid line is a numerical fit to the data. The error bars are smaller than the size 
of the symbols. 

Fig. 11. Isothermal compressibility of liquid ^He as a function of the density. The same 
notation as in Fig. 10. 
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E/N 


T/N 


V/N 


VMC-V-ji 


-5.683±0.014 


15.119±0.005 


-20.802±0.009 


VMC-'0j2 


-5.881±0.005 


15.248±0.004 


-21.129±0.007 


YMC-tpjT 


-6.617±0.007 


14.552±0.030 


-21.169±0.018 




7 11 c^_i_n ni n 
- < .iiOitU.UiU 


1 A c;cn_i_n non 
i4.0oy±U.UzU 


01 vn/i_i_n non 
-zi. ; U4±U.UzU 


QDMC-V'J2 


-7.121zb0.010 


14.576±0.025 


-21.697±0.023 


QDMC-V' JT 


-7.125±0.005 


14.417±0.030 


-21.542±0.020 


GFMC 


-7.120±0.024(") 


14.47±0.09('') 


-21.59±0.09('') 



Table I 
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A.ziz 






Aziz II 






p{a ^) 




E/N 




/ AT 

T/N 


V/N 


h/N 


T/N 


V/N 


E/N 


0.328 


-6 


988±0 


013 


12.107±0.018 


-19.095±0.013 


-7.150±0.010 


12.152±0.032 


-19.302±0.030 






(-7 


034±0 


037) 














0.365 


-7 


121±0 


010 


14.576±0.025 


-21.697±0.023 


-7.267±0.013 


14.622±0.027 


-21.889±0.024 


-7.17 




(-7 


120±0 


024) 














0.401 


-6 


892±0 


013 


17.262±0.030 


-24.154±0.027 


-7.150±0.016 


17.302±0.038 


-24.452±0.035 


-7.03 






Qq4_|_r) 


048 "1 














0.424 


-6 


696±0 


024 


19.152±0.042 


-25.848±0.035 


-6.877±0.022 


19.218±0.037 


-26.095±0.030 


-6.77 


0.438 


-6 


422±0 


020 


20.447±0.036 


-26.869±0.030 


-6.660±0.017 


20.398±0.034 


-27.058±0.030 


-6.55 




(-6 


564±0 


058) 














0.490 


-5 


OlOzbO 


025 


25.402±0.047 


-30.412±0.040 


-5.222±0.025 


25.404±0.050 


-30.626±0.043 






(-5 


175±0 


101) 















Table II 
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P(atm) 


K(atm ^) 


c(m/s) 


Aziz 


0.878±0.073 


0.01199±0.00004 


241.53±0.44 


Aziz II 


-0.019±0.075 


0.01241±0.00004 


237.40±0.46 


Exp 


0. 


0.0124 


237.2 



Table III 
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Aziz II + AT 


Aziz II + BM 


Aziz II -1- MBM 


Exp 




E/N 


V/N 


E/N 


V/N 


E/N 


V/N 


E/N 


0.328 


-7.045 


-19.197 


-7.122 


-19.274 


-7.071 


-19.223 




0.365 


-7.127 


-21.749 


-7.249 


-21.871 


-7.168 


-21.790 


-7.17 


0.401 


-6.971 


-24.273 


-7.141 


-24.443 


-7.027 


-24.329 


-7.03 


0.424 


-6.668 


-25.886 


-6.886 


-26.104 


-6.741 


-25.959 


-6.77 


0.438 


-6.435 


-26.833 


-6.675 


-27.073 


-6.515 


-26.913 


-6.55 


0.490 


-4.913 


-30.317 


-5.323 


-30.727 


-5.050 


-30.454 
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